Simulated Examples

Constant Mean (Ordinary Kriging)

First, we simulate data from a GRF with mean 0 and an exponential covariance structure with parameters \(\sigma^2 = 5\) and \(\phi = 30\). We simulate observations on a sample of size 600 from the set \(D = \{1, \ldots, 100\} \times \{1, \ldots, 100\} \subset \mathbb{R}^2\).

library("tidyverse"); theme_set(theme_bw())
library("patchwork")

set.seed(2)

grid <- expand.grid(1:100, 1:100) 

dom <- grid %>% 
  tibble() %>% 
  rename(x = Var1, y = Var2) %>% 
  slice_sample(n = 600)

data <- geoR::grf(grid = dom, cov.model = "exp", cov.pars = c(5, 30), nugget = 1)$data %>%
  {tibble(z = .)}

df <- bind_cols(dom, data)

Below, we have visualized the simulated data and plotted the empirical variogram.

ggplot(df, aes(x, y, color = z, size = z)) +
  geom_point() +
  scale_color_viridis_c() +
  guides(color = guide_legend(), size = guide_legend())

vgm <- geoR::variog(data = df[,3], coords = df[,1:2])
plot(vgm)

We can estimate parameters of the variogram and perform ordinary kriging to interpolate between the original data. Note, we have to provide initial values for the parameters to be estimated as well as specify the exponential covariance model.

vgm_fit <- geoR::variofit(vgm, ini.cov.pars = c(6, 40), cov.model = "exp")

krige <- geoR::krige.conv(data = df$z, coords = df[,1:2], locations = grid, 
                    krige = geoR::krige.control(obj.model = vgm_fit))

df_krige <- bind_cols(grid[,1], grid[,2], krige$predict) %>% 
  rename(x = `...1`, y = `...2`, z = `...3`)

ggplot(df, aes(x, y, size = z)) +
  geom_raster(data = df_krige, aes(x, y, fill = z)) +
  geom_point() +
  coord_cartesian(xlim = c(1, 100), ylim = c(1, 100), expand = FALSE) +
  scale_fill_viridis_c(name = "prediction") +
  guides(size = guide_legend(reverse = TRUE))

Now, we load in the functions which extend ggplot2:

source(here::here("ggplot_work/Spatial_Modeling/R/geom-krige.R"))
source(here::here("ggplot_work/Spatial_Modeling/R/stat-krige.R"))
source(here::here("ggplot_work/Spatial_Modeling/R/geom-krige-contour.R"))
source(here::here("ggplot_work/Spatial_Modeling/R/stat-krige-contour.R"))

Using these functions, we can create the same visualizations much easier. Below, see that the kriging is handled by geom_krige – we only need supply the model specification and initial values for the parameters. Recall, these can be determined by visualizing the empirical variogram as before.

ggplot(df, aes(x, y)) +
  geom_krige(aes(z = z), inits = c(6, 40), model = "exp") +
  geom_point(aes(size = z)) +
  coord_cartesian(xlim = c(1, 100), ylim = c(1, 100), expand = FALSE) +
  scale_fill_viridis_c() +
  guides(size = guide_legend(reverse = TRUE))

We can adjust the resolution of the kriging as shown below:

p_low <- ggplot(df, aes(x, y, z = z)) +
  geom_krige(nx = 50, ny = 50, inits = c(6, 40), model = "exp") +
  scale_fill_viridis_c() +
  coord_cartesian(xlim = c(1, 100), ylim = c(1, 100), expand = FALSE) 

p_high <- ggplot(df, aes(x, y, z = z)) +
  geom_krige(nx = 200, ny = 200, inits = c(6, 40), model = "exp") +
  scale_fill_viridis_c() +
  coord_cartesian(xlim = c(1, 100), ylim = c(1, 100), expand = FALSE)

p_low / p_high

We can also visualize the kriging output via contours with geom_krige_contour and geom_krige_contour_filled (which function similarly to geom_contour and geom_contour_filled):

p_contour <- ggplot(df, aes(x, y, z = z)) +
  geom_krige_contour(inits = c(6, 40), model = "exp") +
  coord_cartesian(xlim = c(1, 100), ylim = c(1, 100), expand = FALSE)

p_contour_filled <- ggplot(df, aes(x, y, z = z)) +
  geom_krige_contour_filled(inits = c(6, 40), model = "exp") +
  scale_fill_viridis_d() +
  coord_cartesian(xlim = c(1, 100), ylim = c(1, 100), expand = FALSE)

p_contour | p_contour_filled

We have also allowed for control of the resolution of the kriging here:

p_contour_lower <- ggplot(df, aes(x, y, z = z)) +
  geom_krige_contour(inits = c(6, 40), model = "exp", nx = 30, ny = 30, bins = 5) +
  coord_cartesian(xlim = c(1, 100), ylim = c(1, 100), expand = FALSE)

p_contour_filled_lower <- ggplot(df, aes(x, y, z = z)) +
  geom_krige_contour_filled(inits = c(6, 40), model = "exp", nx = 30, ny = 30, bins = 5) +
  scale_fill_viridis_d() +
  coord_cartesian(xlim = c(1, 100), ylim = c(1, 100), expand = FALSE)

p_contour_lower | p_contour_filled_lower

Non-Constant Mean (Universal Kriging)

The extensions also allow for the mean of the GRF to depend on the coordinates linearly. Below, we simulate from such a process using the same set \(D\) as before. However, this time we are instead simulating from a process with Gaussian covariance and different covariance parameters.

set.seed(1)
data <- geoR::grf(grid = dom, cov.model = "gau", cov.pars = c(10, 20), nugget = .5)$data %>%
  {tibble(z = .)}

df <- bind_cols(dom, data) %>% 
  mutate(z = z + 3/20 * x - 1/20 * y)

Below, we can see the data is spatially corellated in addition to the obvious linear trend:

ggplot(df, aes(x, y, color = z, size = z)) +
  geom_point() +
  scale_color_viridis_c() +
  guides(color = guide_legend(), size = guide_legend())

In order to allow for this nonconstant mean, we perform Universal Kriging. This is available in the extension functions via the argument method = "UK":

geoR::likfit(data = df$z, coords = df[,1:2], trend = "1st", ini.cov.pars = c(1, 1), cov.model = "gau")
ggplot(df, aes(x, y)) +
  geom_krige(aes(z = z), method = "UK", inits = c(1, 1), model = "gau") +
  geom_point(aes(size = z)) +
  coord_cartesian(xlim = c(1, 100), ylim = c(1, 100), expand = FALSE) +
  scale_fill_viridis_c() +
  guides(size = guide_legend(reverse = TRUE))

Below, we perform ordinary kriging and ignore the spatial structure of the mean. See that the resulting interpolation does not fit the data as well, missing out on local variation.

ggplot(df, aes(x, y)) +
  geom_krige(aes(z = z), method = "OK", inits = c(8, 10), model = "gau") +
  geom_point(aes(size = z)) +
  coord_cartesian(xlim = c(1, 100), ylim = c(1, 100), expand = FALSE) +
  scale_fill_viridis_c() +
  guides(size = guide_legend(reverse = TRUE))

Below, we show each of the different methods of visualizing universal kriging:

p_heatmap_UK <- ggplot(df, aes(x, y)) +
  geom_krige(aes(z = z), method = "UK", nx = 200, ny = 200, inits = c(1, 1), model = "gau") +
  coord_cartesian(xlim = c(1, 100), ylim = c(1, 100), expand = FALSE) +
  scale_fill_viridis_c() 

p_contours_UK <- ggplot(df, aes(x, y)) +
  geom_krige_contour(aes(z = z), method = "UK", nx = 30, ny = 30, inits = c(1, 1), model = "gau") +
  coord_cartesian(xlim = c(1, 100), ylim = c(1, 100), expand = FALSE) 

p_filled_contours_UK <- ggplot(df, aes(x, y)) +
  geom_krige_contour_filled(aes(z = z), method = "UK", nx = 30, ny = 30, inits = c(1, 1), model = "gau") +
  coord_cartesian(xlim = c(1, 100), ylim = c(1, 100), expand = FALSE) +
  scale_fill_viridis_d()

p_heatmap_UK /
p_contours_UK /
p_filled_contours_UK

Note, this relies on the geoR function likfit which requires good initial values. Sometimes, trial and error is required – the function recommends trying numerous initial values. Interestingly, the initial values of \((1, 1)\) seem to perform well, regardless of the true values of the covariance parameters.

Also, it may be possible to specify the processes mean’s relationship with the coordinates via a formula. This would be a good option, it could allow for polynomial relationships.

Latitude Longitude Data and ggmap

These extensions work well with the ggmap function. Below, we simulate data in the Waco area:

set.seed(3)

df <- tibble(
  lon = runif(250, -97.3, -97),
  lat = runif(250, 31.4, 31.7),
  z = geoR::grf(grid = cbind(lon, lat), cov.model = "exp", cov.pars = c(3, .1), nugget = 1)$data
)

Below, we set up ggmap and visualize our sample points:

library("ggmap")

theme_update(
  axis.title = element_blank(),
  axis.text = element_blank(),
  axis.ticks = element_blank() 
)

waco_map <- get_map(location = c(-97.3, 31.4, -97, 31.7), color = "bw")

ggmap(waco_map) +
  geom_point(data = df, aes(lon, lat, color = z, size = z)) +
  guides(color = guide_legend(), size = guide_legend()) +
  scale_color_viridis_c()

It is simple to include kriging layers on top of ggmap raster images, we need only adjust alpha:

p_heatmap_dots <- ggmap(waco_map) +
  geom_krige(data = df, aes(lon, lat, z = z), alpha = .4, inits = c(5, .5), model = "exp") +
  geom_point(data = df, aes(lon, lat, size = z)) +
  scale_fill_viridis_c()

p_filled_contours_dots <- ggmap(waco_map) +
  geom_krige_contour_filled(data = df, aes(lon, lat, z = z), alpha = .4, bins = 7, 
                            inits = c(5, .5), model = "exp") +
  geom_point(data = df, aes(lon, lat, size = z)) +
  scale_fill_viridis_d()

p_heatmap_dots /
p_filled_contours_dots

And now, without the plotted observations:

p_heatmap <- ggmap(waco_map) +
  geom_krige(data = df, aes(lon, lat, z = z), alpha = .4, inits = c(5, .5), model = "exp") +
  scale_fill_viridis_c()

p_filled_contours <- ggmap(waco_map) +
  geom_krige_contour_filled(data = df, aes(lon, lat, z = z), alpha = .4, bins = 7, 
                            inits = c(5, .5), model = "exp") +
  scale_fill_viridis_d()

p_contours <- ggmap(waco_map) +
  geom_krige_contour(data = df, aes(lon, lat, z = z), alpha = .4, bins = 7, 
                     inits = c(5, .5), model = "exp")

p_heatmap /
p_filled_contours /
p_contours